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ABSTRACT 

We present a mass model for the Milky Way, which is fit to observations of the Sagit- 
tarius stream together with constraints derived from a wide range of photometric 
and kinematic data. The model comprises a Sersic bulge, an exponential disk, and an 
Einasto halo. Our Bayesian analysis is accomplished using an affine-invariant Markov 
chain Monte Carlo algorithm. We find that the best-fit dark matter halo is triaxial with 
axis ratios of 3.3 ±0.7 and 2.7 ±0.4 and with the short axis approximately aligned with 
the Sun-Galactic centre line. Our results are consistent with those presented in Law 
and Majewski (2010). Such a strongly aspherical halo is disfavoured by the standard 
cold dark matter scenario for structure formation. 
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1 INTRODUCTION 

It is a prediction of the hierarchical clustering scenario, 
borne out by N-body simulations, that dark matter halos 
are triaxial (see, for example, Frenk (1988), Franx, Illing- 



worth, fe de Zeeuw (1991), Warren et al. (1992), Jing fc 



Suto (2002), and Allgood et al. (2006)). To be sure, baryons 



will alter the shapes of halos, presumably making them more 
spherical (Dubinski 1994; Kazantzidis et al. 2004; Debattista 



et al.|[2008 ). Nevertheless, halo triaxiality is expected to be 
a feature of our Universe if structure formation proceeds ac- 
cording to the standard model. Unfortunately, the shapes 
of dark matter halos are extremely difficult to determine by 
direct observations. 

Tidal debris from the Milky Way's satellite galaxies pro- 
vide a potentially powerful probe of the Galactic gravita- 
tional potential with the tidal stream from the Sagittarius 
dwarf spheroidal galaxy (Sgr dSph) the most prominent ex- 
ample. In the years following the discovery of the Sgr dSph 



( jlbata et al.|1997| ) numerous surveys have produced an all- 
sky panorama of the Sgr system (see [Majewski et aL] (|2003) 
and references therein). The observed positions and veloc- 
ities of stars in the stream are compared with theoretical 
predictions for a range of triaxial halo models. The under- 
lying assumption is that the kinematics of the stream is 
determined primarily by the Galactic potential, and only to 
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a lesser extent, by the structure of the progenitor. Examples 
of studies that attempted to use the Sgr system in this way 



include Ibata et al. (2001), Helmi (2004), Martmez-Delgado 


et al. (2004 


, Johnston, Law, &; Majewski (2005), and Fell- 


hauer et al. 


(2006 


)• 



Majewsk i et al.| ( |2003| presented an all-sky map of the 



Sgr system in M-giant stars in which candidate Sgr stars 
were selected from the 2MASS catalog using their color, 
magnitude, and angular position. The result was a kinematic 



snapshot of both leading and trailing streams. Bel okurov et| 
|aL| (|2006[) provided a more refined map of the stream us- 



ing the SDSS. Recently Law, Majewski, & Johnston (2009) 



and |Law fe Majewski| ( |2010| (hereafter LMJ09 and LM10, 
respectively) used both of these data sets to argue that the 
Milky Way's halo is triaxial. LMJ09 compared single par- 
ticle orbits to the projected positions of the SDSS fields 
and radial velocities of the M giants. Their best-fit model 
has isopotential axis ratios of = 1.5 and q z ^ = 1.25 
with the halo's intermediate axis pointing toward the North 
Galactic pole and its minor axis roughly aligned with the 
line that connects the Sun and the Galactic centre (GC). 
LM10 carried out a similar analysis, but modeled Sagittar- 
ius as an N-body system rather than a single particle. Their 
favored halo model is nearly oblate with ~ q z ,<s> — 1.4 
where again the minor axis of the halo lies in the disk plane 
nearly along the Sun-GC line. 

The LMJ09 and LM10 results present a challenge for the 
standard cold dark matter model of structure formation. In 
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general, a halo's gravitational potential will be more spher- 
ical than its mass distribution; the isodensity and isopoten- 
tial axis ratios of an axisymmetric system (q p and q<s> re- 
spectively) roughly satisfy the relation q p ~ 3 (q<z> — 1) + 1, 



at least in the outer parts of the halo (Binney & Tremaine 



2008). Thus, the mass distribution that corresponds to the 



LM10 potential is approximately oblate with axis ratios 
~ 2.2. Since the formation of the Galactic disk is expected 
to make the halo more symmetric in the disk plane the axis 
ratios of the proto-galaxy are likely even larger. Allgood et 
al. (2006) presented a comprehensive analysis of dark halo 
shapes based on high-resolution dissipationless simulations. 
They found that the mean major-to-minor isodensity axis 
ratio for Galaxy-sized halos was q p ~ 1.7 with a dispersion 
of a q ~ 0.3. Moreover their halos, and those of earlier stud- 
ies, tended more toward prolate-triaxial than oblate-triaxial 
(that is triaxial with one axis significantly smaller than the 
other two). The effect baryons have on a halo's shape was 



studied by |Dubinski (1994); Kazantzidis et al. (2004); De- 



battista et al. (2008) and others. In general baryons tend to 
make the halos more spherical, particularly in the innermost 
regions. The axis ratios of the proto-galaxy should therefore 
be even larger than those observed today. The implication 
is that the Galactic halo favoured by the LMJ09 and LM10 
analyses is an outlier within the context of the standard 
cosmological paradigm. 

A secondary issue concerns the shapes of the isoden- 
sity contours themselves. LMJ09 and LM10 assume triaxial 
isopotential contours, which imply peanut-shaped isoden- 
sity contours (see, for example, |Binney &; Tremaine| ( |2008| ) . 
By contrast, triaxiality in simulations is almost always mea- 
sured in terms of the density and peanut-shaped contours 
are rarely seen. 

These considerations motivated us to reexamine con- 
straints on the dark halo from the Sgr stream. We carry 
out the modeling exercise within the Bayesian statistical 
framework and are therefore able to calculate the prob- 
ability distribution function (PDF) for the structural pa- 
rameters of the halo and quantify the model uncertainties. 
Previous analyses focused on the shape parameters of the 
halo (the axis ratios and the Euler angles that define the 
halo's orientation). The general practice has been to fix the 
structural parameters of the bulge, disk, and spherically- 
averaged halo and allow the halo shape parameters to vary 
while trying to fit kinematic data for stars along the Sgr 
stream. Typically, only the most rudimentary constraints 
on the Galaxy's structure are included. For example, in the 



LMJ09 and LM10 analyses (see, also | Johnston et al.| (|2005)) 
the model is designed so that the circular speed at the posi- 
tion of the Sun is v c — 220kms _1 . As discussed below, the 
disk and bulge in this model appear to be too large (by fac- 
tors of 2-3) to be consistent with a number of observational 
constraints. 

In this paper, we consider a general disk-bulge-halo 
model for the Milky Way and allow the key structural pa- 
rameters for all three components to vary simultaneously. 
The likelihood function includes not only the stream con- 



straints, but observational constraints from the line-of-sight 
velocity dispersion (LOSVD) and surface brightness profiles 
for the bulge, the vertical force and surface density in the 
solar neighborhood, the Oort constants, the circular speed 
curve, and the mass at large Galactocentric radii. 

At first glance, our parameter space has expanded to 
an unwieldy level. (Our model has twenty-eight free pa- 
rameters!). However, by using a Markov chain Monte Carlo 
(MCMC) algorithm we are able to efficiently estimate the 
PDF for the full multi-dimensional parameter space. In this 
paper, we employ the affine-invariant 'Stretch-Move' ensem- 



ble sampler of Goodman & Weare (2010) (see also | Foreman 
|Mackey et al. (2012)). Since our MCMC analysis requires a 



large number of "likelihood calls" we are precluded from 
modeling Sgr as a full N-body system. We therefore follow 
LMJ09 and model Sgr and the associated tidal stream as a 
single particle that orbits in the fixed gravitational potential 
of the Milky Way. The stream likelihood function is calcu- 
lated by comparing the phase space coordinates from the M 
giant and SDSS observations (namely, the angular positions 
and line-of-sight velocities) with points along model orbits. 

In Section 2 we use Bayesian inference and our MCMC 
algorithm to reanalyze the M giant and SDSS data within 
the context of the Galactic model from LMJ09 and LM10. 
We also consider an alternative model in which the isoden- 
sity contours (rather than isopotential contours) of the halo 
are triaxial. In Section 3 we introduce a more general model 
for the density-potential pair of the Galaxy as well as the 
observational constraints that enter our analysis. We present 
our main results in Section 4. Our analysis points to a sim- 
ilarly triaxial model as the one found in LM10 though the 
disk-halo decomposition of the mass distribution is quite dif- 
ferent. We conclude in Section 5 with a summary of our key 
results. 



2 HALO TRIAXIALITY FROM THE 

SAGITTARIUS STREAM VIA BAYESIAN 
INFERENCE 

We begin this section with brief summaries of the Galac- 
tic model used in LMJ09 and LM10, the observational con- 
straints from the Sgr stream, and the affine-invariant en- 
semble sampler deployed in our analysis. We then present 
results from two separate MCMC runs. The first closely fol- 
lows LMJ09 in that the isopotential contours are assumed to 
be triaxial ellipsoids and certain model parameters, namely 
the present-day distance to the Sgr dSph and the thickness 
and intrinsic dispersion of the stream are kept fixed. In the 
second run, the isodensity contours are assumed to be triax- 
ial and the aforementioned parameters are allowed to vary. 



2.1 Galactic model 

The Galactic model used by LMJ09 and LM10 comprises a 
Hernquist bulge (|Hernquist 1990), a Miyamoto-Nagai disk 



(Mi yamoto &; Nagai||1975 ), and a logarithmic halo (see, for 
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example, Binney & Tremaine (2008)). The contributions to 



the gravitational potential for each of these components are 
given by 

GM b 



r + c 



GM d 



^R 2 + (a + Vz 2 + b 2 ) 2 
Q h (rt) = a 2 h \n(r 2 + d 2 ) . 



(1) 
(2) 
(3) 



where (x, y, z) are the usual right-handed Cartesian coordi- 
nates with the z-axis oriented along the symmetry axis of the 
disk and the Sun located on the x-axis, R — (x 2 + y 2 ) 1 ^ 2 , 



2\V2 



We also define a triaxial coordinate 



and r = [R 2 + 
system for the halo with coordinates r t = (x*, yt, zt). Fol- 
lowing LMJ09, we fix the z^-axis to be along the symme- 
try axis of the disk and write rt = RA$r where A$ = 
diag(l,A- 1 ,B- 1 ) and 



R 



cost/ 
— sin^ 





(4) 



In this section, we follow LMJ09 (see also Johnston et al 



(1999) and Law, Johnston fe Majewski (2010)) and fix the 



structural parameters of the disk, bulge, and halo potentials 
to the following values: M b = 3.43 x 1O 1O M , c = 0.7 kpc, 
M d 1.0 x 10 n M o , a 6.5 kpc, b 0.26 kpc, and d 
12kpc. The halo scale velocity, an, is adjusted so that the 
local circular speed is 220 km s _1 . We refer to this model as 
the LJM model. 



2.2 Observational constraints from the Sgr dSph 
and stream 

The Sgr system roughly lies in a single plane, which has been 
identified as the orbital plane of the Sgr dSph. To a good 
approximation, this plane contains the Sun and is defined 
by the orbital pole b p ) = (273.8°, -14.5°) (LM10 and 
| Johnston et al.| ([2005)). The Sgr dSph itself is located at 
(Zsg r, 6 Sgr ) = (5. 6°, -14.2°). 

|Ibata et al.| p997) found that the stars observed in the 
central regions of the Sgr dSph have a mean radial velocity 
of 171 =b 1 km s _1 , which is nearly the same as the radial ve- 
locity of M54, a globular cluster usually identified with the 
center of the Sgr dSph. Note that this velocity has been con- 
verted to the Galactic Standard of Rest (GSR), a reference 
frame centered on the Sun but at rest relative to the Galac- 
tic center. The conversion from heliocentric velocities to the 
GSR depends on the circular speed at the position of the 
Sun and the Sun's peculiar velocity (£/©, V©, W©). In |Ibata| 



et al. (1997) as well as LMJ09 and LM10, the circular speed 



of the Sun is assumed to be 220 km s _1 . We return to this as- 
sumption in Section 3. The heliocentric proper motion of the 
Sgr dSph has been measured by |Pryor, Piatek, Olsweski 
( 2007J) using archival Hubble Space Telescope (HST) data as 



(fi h fi b ) = (-2.615 ± 0.22, 1.87 ± 0.19) masyr" 1 . The helio- 
centric distance to the Sgr dSph i^sgr has been estimated 
by a various methods and is found to lie in the range from 
22 - 28.4 kpc (see Table 2 of |Kunder fc Chaboyer| | |2009) ). 



2.3 Likelihood function and posterior probability 
distribution 

For simplicity we use the so-called Sagittarius spherical co- 



ordinate system (d, A, f3) defined by Majewski et al. (2003) 
where d is the heliocentric radial coordinate and A and f3 are 
angular coordinates. The Sgr dSph is located at A = /3 = 0° 
while the stream approximately follows the f3 — 0° great 
circle with A increasing along the trailing portion of the 
stream and f3 increasing toward the orbital pole. Follow- 
ing LMJ09, we assume that the Sgr dSph can be modeled as 
a point mass and that the leading and trailing portions of 
the stream follow, respectively, the future and past segments 
of its orbit. The model orbit is then compared with the ra- 
dial velocity and angular position measurements of M giants 



from Majewski et al. (2003) and the angular position mea- 



surements from the SDSS study of Belokurov et al. (2006). 



This method ignores the internal structure and dynamical 
evolution of the progenitor and stream. N-body simulations 
allow one to explore these effects but their use is unfeasi- 
ble in the present analysis where 500K likelihood calls are 



required. (See however Varghese et al. (2011) who discuss 
a promising approach in which a family of kinematic orbits 
are used to model the "thick" tidal stream.) 

For the "i'th" data point, we determine the position 
along the orbit with the same A^ and, in so doing, arrive 
at model predictions for the line-of-sight velocity vm (A*) 
and the angular position /3m (A«). The likelihood function, 
or equivalently, the probability of the data given the model 
M, is then 



'(^w = n 7^172 — exp 

f = i (2tt) 1 (J^ % 



(Pm (A») - Pi) 2 



N v 

x n^i75 



exp 



2 °h 



(VM (Ai) - Vi) 



2a, 



(5) 



The first term on the right-hand side involves a product over 



all M giants in the Majewski et al 



(2003) data set as well 



as the SDSS Sgr stream fields from |Belokurov et al.| (|2006) 
whereas the second term involves only the M giants. The 
parameters ap are meant to account for both observational 
uncertainties and the intrinsic thickness of the Sgr stream. 
Likewise, the <j v are meant to account for observational un- 
certainties in the line-of-sight velocity and the intrinsic dis- 
persion. 

Our aim is to calculate the posterior probability func- 
tion (PDF) of the model, which is given by Baye's theorem, 



p(M\D,I) 



p(M\I)p(D\M) 



(6) 



p(D\I) 

where I represents prior information and P(M\I) is the prior 



© 0000 RAS, MNRAS 000, 000-000 



4 N. Deg & L. Widrow 



probability on the model M. The term p(D\I), often referred 
to as the evidence, is essentially a normalization factor and 
does not enter our calculations. 

In general, M is specified by Np parameters and there- 
fore P(M\D, I) is an A/p-dimensional function. In order 
to efficiently map this multi-dimensional function we use 
the Stretch-Mov e (SM) MCMC algorithm from [Goodman 
& Weare|(|2010|). All MCMC algorithms require a proposal 



distribution or sampler to generate the chain of points in 
parameter space. If the different model parameters have dif- 
ferent scales or are strongly correlated, then it can be ex- 
tremely difficult to choose an effective sampler. And a poorly 
chosen sampler leads to poor convergence of the chain. The 
SM-MCMC algorithm employs an ensemble of Nw "walk- 
ers" where Nw > Np . The proposal distribution for a given 
walker is determined by the other walkers (rather than by 
the user). By design, the algorithm is invariant under linear 
re-parameterizations of the model parameters (i.e., afflne 
invariant) and the difficulties of choosing a sampler are 
avoided. 



2.4 LMJ09 revisited 



We next reexamine the Majewski et al. ( 2003]) and |Belokurov| 
et al.| ( |2006| data sets within the context of the Galactic 
model described above and model assumptions that are sim- 
ilar to those used in LMJ09. The free parameters are the axis 
ratios A<$> and the angular position of the x^-axis in the 
xt — yt plane and the two proper motion components of the 
Sgr dSph. Following LMJ09, we fix Ds g r and line-of-sight ve- 
locity for the Sgr dSph to be 28kpc and 171 km s" 1 (GSR), 
respectively. We also follow LMJ09 and set ap^ = 1.9° and 
a Vj i — 12 km s _1 . 

We assume uniform priors in log (A) and \og(B) since 
the axis ratios 1 : 0.5 and 1 : 2 are effectively equivalent. To 
avoid degeneracies, we restrict A to be greater than unity. 
The prior for is assumed to be uniform between 0° and 
180°. Finally, we assume that the prior probability on the 
components of the proper motion are uniform across a range 
that includes both the best-fit values from LM10 and values 



observed by |Pryor et al. ( 2007 ) . 

We fit the stream using our SM-MCMC algorithm with 
100 walkers and 5000 steps. Here and throughout this work 
we discard the first 1000 steps, i.e., the burn-in segment of 
the chain. The marginalized PDFs for B®, and 6 in 
this run (referred to as the triaxial potential or TP run) are 
shown in Figure 1 while a summary is provided in Table 
1. Interestingly, the most likely shape from the Bayesian 
analysis is closer to the results of LM10 than the results 
from LMJ09. The implication is that the parameter search 
algorithm is at least as important as whether the Sgr dSph 
is modeled as an N-body system or a point particle. 

2.5 Triaxial density 

By design, isopotential surfaces in the halo model described 
by Equation [3] are triaxial whereas the isodensity surfaces are 



not. Indeed, when A and B differ significantly from 1, as in 
the best-fit models of LMJ09 and LM10, isodensity surfaces 
become peanut-shaped (see Fig. 2.9 of |Binney &; Tremaine| 
(12008}). In simulations, the shape of dark matter halos is 
inferred by measuring the axis ratios of isodensity surfaces 
assuming that these are the surfaces that are triaxial. With 
these issues in mind, we consider an alternative model in 
which the halo density is given by 



Ph (R,z)= Ph (r 2 )^ 



2 r 2 + 3d 2 



2nG 1 



+ <py 



(?) 



where r t = RA p r with A p = diag (l, A" 1 , B" 1 ) . Note that 
in the spherical limit, r = r t and Equations [7] and [3] describe 
equivalent models. 

The force and potential are calculated from Equation [7] 
via the homeoid theorem (|Binney fc Tremaine 2008). The 



potential for any triaxial density is given by 

ip(oo) — ip(m) 



ai 



dr- 



y/(T + al)(T + 0*)(T + d 

where the ai are the axis ratios, 



2 



a? +r 



is similar to the square of the ellipsoidal radius, and 



dm 2 p(m 2 ) 



(8) 



(9) 



(10) 



1, CL2 



is an auxiliary function. In our models a\ 
a 3 B p . 

We consider a Galactic model in which Equation [7] 
rather than Equation [3] describes the halo. We also allow 
Dsgr to vary within the range of measured values compiled 
by |Kunder fe Chaboyerj fl2009| ). To be precise, we assume a 
uniform prior for Ds gT between 22kpc and 28.5 kpc. |Law et| 
al. (2010) found that the effect of varying Ds g r was degen- 



erate with other parameters, particularly the distance R Q of 
the Sun to the Galactic center. We therefore include R as 
a model parameter and assume a Gaussian prior centered 



on 8.2 kpc with 1-a width of 0.4 kpc (Bovy, Hogg, & Rix 



[2009} . We also include a prior probability for the proper 



(2007) 



motion based on the HST analysis of |Pryor et al 
Finally, we replace ap^ and <j v ,i in Equation [5] with the free 
parameters epj and e v , where j — MG for the [ Majewski 



et al. 



d2003[ ) M giants and SDSS for the |Belokurov et al. 



(2006) SDSS fields. That is, we model the stream thickness 
in both angular position and line-of-sight velocity. To sum- 
marize, our model now includes five additional parameters, 
£*s g r, Ro, £mg, csdss, and e v . 

The marginalized PDFs for A p , B p , and are shown in 
Figure 2 (referred to as the triaxial density or TD run) . We 
also compare the best-fit values with those from the previous 
analysis and with the results from LMJ09 and LMJ10. To 

3 (A$ — 1) + 1 from Binney 



do so, we use the relation A p 



& Tremaine (120081). Once again, we find that the Sgr stream 



favours a strongly oblate-triaxial halo with isodensity axis 
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A<$> 


-£><!> 


A P 


Bp 


e 


LMJ09 


1.5 


1.25 


(2.5) 


(1.75) 


0° 


LM10 


1.4 


1.4 


(2.2) 


(2.2) 


-7° 


TP 


1.49 


1.34 


(2.48) 


(2.03) 


-9° 


TD 


(1.49) 


(1.38) 


2.48 


2.15 


-11° 



Table 1. Axis ratios and angle 6 for the LMJ09, LM10 analyses and TP and TD runs as 
described in Section 2. Parenthesis denote axis ratios that are calculated using the expression 
A p = 3 (A$ — 1) + 1 or its inverse. 



ratios that are greater than 2 and with the symmetry axis 
closely aligned with the Sun-GC line. The uncertainty in the 
axis ratios is greater, by about a factor of three, for this run. 

In Figure 3, we show the marginalized PDF in the 
Dsgr — Ro plane. The constraint on R Q comes almost entirely 
from the prior, while our fit to the stream data favours the 
upper range of observed values for the distance to the Sgr 
dSph. 

In Figure 4, we show the marginalized PDF for the com- 
ponents of the proper motion, pi and pb- Evidently, the 
PDFs are inconsistent with measurements by |Pryor et al.| 
(2007) at about the 2cr-level but are consistent with the 



favoured values obtained in LMJ09 and LM10. Tension be- 
tween the model and observed values for the proper motion 
may indicate a problem with the model assumptions (e.g., 
interpretation of the streams as direct tracers of the orbit) 
though clearly a more refined measurement of pi and pb is 
required. 



3 TOWARD A MORE REALISTIC MILKY 
WAY MODEL 

3.1 An updated Galactic model 

The model described by Equations 1-3 is attractive in large 
part because of its simplicity since the potential, force, and 
density can all be expressed in terms of analytic functions. 
However, this model does not represent our current un- 
derstanding of the azimuthally-averaged structure of spi- 
ral galaxies. For example, Andredakis, Peleti er7 & BalceHs] 
( |1995| showed that bulges of spiral galaxies could be best 
fit by the Sersic law. Furthermore, since the classic paper by 
Freeman (1970), it has become standard practice to model 
the luminosity profile of disks as an exponential. Finally, 



the spherically-averaged density profile of dark matter ha- 



los is described extremely well by the Einasto profile (Mer- 



ritt et aL]|2006| . Therefor we consider a Milky Way model 
with these three components. We assume that the bulge and 
disk have constant (but different) mass-to-light ratios. The 
density that produces a Sersic surface brightness profile is 
flPrugniel fe Simien|[l997| jTerzic fe Graham|2005t 



Pb(r) p b0 



(£)"*■ 



(ID 



where p = 1 — 0.6097/n + 0.05563/n 2 , n is the Sersic index 
and b — b(n) is chosen so that the radius R e encloses half 



the total projected mass. For simplicity, we parameterize the 
bulge by the scale velocity 

a b = (4irnb n(p - 2) r (n (2 - p)) R 2 ePb0 ) 1/2 (12) 

rather than pbo. 

We assume that the disk density is exponential in the 
radial direction (the Freeman Law) and has a sech 2 structure 
in the vertical direction: 



p d (R, z) 



M d 



4nR 2 d z d 



~ R/Rd sech 2 (z/z d ) 



(13) 



where M d , Rd, and Zd are the mass, radial scale length, 
and vertical scale height for the disk, respectively. The disk 
potential is calculated using the technique of Kuijk en~"fe] 
Dubinski (1995). An analytic 'fake' disk density-potential 
pair, (p/d, &fd) is constructed so that pd — pfd + pr and 
&d = <&fd + where (p r , <E> r ) is the density-potential pair 
of the residual. The fake disk is designed to account for the 
high-order moments of the disk potential. The full poten- 
tial is calculated by solving the Poisson equation using a 
spherical harmonics expansion and iterative scheme. 

The density distribution for our halo model is given by 



Ph(r t ) = poe 



H(rt/r h ) a -l) 



(14) 



where Tt is the triaxial radius, po is a scale density, ru is 
the scale radius, and a controls the logarithmic slope of the 
density profile. The potential for this triaxial density is cal- 
culated via Equation [8] 

To be sure this model represents an approximation 
to the mass distribution and gravitational potential of the 
Galaxy. The model disk and bulge are axisymmetric whereas 
the actual bulge is triaxial and the disk contains a bar. It 
is reasonable to ask whether these non-axisymmetric struc- 
tures could mimic the effects of a triaxial halo. We note that 
at pericenter the Sgr dSph is ~ 20 kpc from the Galactic 
center whereas the bulge triaxiality and the bar are proper- 
ties of the inner few kpc. At pericenter, the disk and bulge 
contribute ~ 25% and 10% respectively to the gravitational 
force. Furthermore, the components of the potential asso- 
ciated with the bar and bulge triaxiality are quadruple (or 
higher order) terms, which fall off significantly faster than 
the leading monopole term. We therefore conclude that the 
triaxiality of the bulge and the presence of a bar will have 
little effect on the orbit of the Sgr dSph and its stream. 
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3.2 Observational constraints 

Our analysis combines kinematic and photometric observa- 
tions that constrain the structure of the Milky Way across 
five orders of magnitude in radius. Apart from the con- 
straints due to the Sgr stream, the observations we con- 
sider are similar to those found in (Dehnen & Binney|1998 



with similar studies (see, for example, Holmberg & Flynn 



Widrow fe Dubinski|2005l |Widrow, Pym, fe Dubinski| 2008 
with some notable exceptions. The starting point is a gen- 
eralization of Equation [5] to 



P (d\s) = n n 



27T<Jij 



(15) 



where the index i runs over the different types of observa- 
tional constraints (e.g. stream angular position and radial 
velocity) while the index j runs over individual data points 
of a particular type of observation. The model prediction is 
M, the observation, O, and is the associated error, a. 

The stream constraints are as before with one proviso: 
Since we are treating the Sun's motion relative to the Galac- 
tic centre as a free parameter, we cannot use the GSR radial 
velocities for the M giants as published in [Majewski et al.| 
(2003) but must use the (observed) heliocentric velocities 
and convert to the GSR using the model values for v c (R ), 
Uq, Vq, and W©. 

For the innermost region of the Galaxy, we use the com- 
pilation of bulge LOSVD measurements found in Tremaine 
et al.| (|2002|), with the restriction that r ^ 4 pc to avoid 



complications from the central black hole. In addition, we 
adjust the dispersion downwards by a factor of 1.07 to ac- 
count for the non-spherical nature of the bulge. We also 
use the surface brightness profile at infrared wavelengths 
from the COBE-DIRBE observations ( |Binney, Gerhard, fe| 
Spergel|1997 ). These observations extend from —40° to 40° 



in Galactic longitude and allow us to probe the structure 
of both the bulge and disk. Of course, our model must now 
include mass-to-light ratios for the bulge and disk. 

We include a number of "local" or solar neighborhood 



constraints in our analysis. Recently Reid et al. (2009) used 



Very Long Baseline radio observations to determine trigono- 
metric parallaxes and proper motions for masers through- 
out the Milky Way's disk. These measurements were then 
used to infer a circular speed at the position of the Sun of 
254 ± 16 km s _1 , a value 15% higher than the one assumed 



in LMJ09 and LM10. The |Reid et al.| ([2009} results were 
scrutinized by |Bovy et al.| ( [2009| ) who carried out a Bayesian 
analysis of the maser data and also incorporated proper mo- 
tion measurement s of Sgr A*. For our analysis, we adopt the 
Bovy et al.|(|2009|) value: v c (R Q ) = 244 ± 13 km s" 1 . 



The shape of the circular speed curve in the solar neigh- 
borhood is described by the Oort constants. We adopt the 
constraints A = 14.8 ± 0.8 km s ^kpc - 1 and B = —12.4 ± 



0.6 km s _1 kpc _1 from |Feast fe Whitelock| ( |l997| ). We use the 
disk surface density measurement of = 49 ± 9 M pc -3 



from Flynn &; Fuchs ( 1994| ) and the vertical force measure^ 
ment of K z (l.l kpc)| (2tt(?) 71 ± 6 M 



©PC 



from 



Kui- 



( [2004} ). 

Observations of the Galactic circular speed curve are 
typically divided into measurements inside and outside the 
solar circle. The inner rotation curve is usually presented 
in terms of the terminal velocity, which is defined as the 
peak velocity along a particular line-of-sight defined by 
the Galactic coordinates 6 = and / where |/| < tv/2. If 
we assume that the Galaxy is axisymmetric, then vterm = 
v c (R)—v c (Ro) sin/. We use data from |Malhotra (1995) with 
the restriction that sin/ ^ 0.3 so as to avoid distortions 
due to the bar. The outer rotation curve requires a more in 
depth discussion. The circular speed v c (R) is related to the 
observed line-of-sight velocity vi sr of a kinematic tracer in 
the disk through the expression 



W(R) = ^v c (R) 



Vc(Ro) = 



Vis 



(16) 



We use observations of HII regions and reflection nebulae 
from|Brand &; Blitz| ([1993 ) and Carbon stars from Demers 
& Battinellil (120071) with the restriction that / < 155° or 



/ ^ 205°, d ^ 1 kpc, and W ^ so as to avoid complications 
due to non-circular motions. The errors in the measurements 
are propagated to errors in W. Additionally, 'noise' param- 
eters for d and vi sr are added in quadrature to the quoted 
errors to account for the potential for unknown systematic 
errors. That is, the velocity error for the j'th data point is 
= (Ty m + where <j v ,m is the measured error (the % sub- 
script has been suppressed) and e u is the noise parameter. 
Similarly, we have a 2 — a 2 + d 2 e 2 d . Note that we allow for 



different noise parameters for the Brand &; Blitz| (|1993) and 
|Demers k Battinell"i| ( |20Q7|) data sets. 



jken & Gilmore (1991). These values are in good agreement 



Finally, we follow |Dehnen Binney ( 1998 ) and use 
M(r < 100 kpc) = (7 ± 2.5) x 10 11 M as a constraint 
on the mass of the Milky Way at large Galacto centric radii. 
This constraint is based on the study of satellite kinematics 
by |Kocha nek (1996), the locally determined 'escape' speed, 
and the timing of the local group. It is also based on mod- 
elling of the Magellanic Clouds and stream by |Lin, Jones, fc| 
Klemola |p995). The large error in M (r < 100 kpc) reflects 
the potential for large systematic errors. In principle, the Sgr 
system should be able to constrain the mass distribution of 
the Milky Way in the 30 - 100 kpc range. 



4 RESULTS 

We now present results based on the constraints described 
in the preceding Section. As before we use the SM-MCMC 
algorithm with 100 walkers and 5000 steps. Our model is 
described by twenty-eight free parameters: nine structural 
parameters for the disk, bulge and halo, mass-to-light ratios 
for the bulge and disk, two axis ratios and an orientation 
angle for the triaxial halo, the distance of the Sun to the 
Galactic center and its motion relative to the local standard 
of rest, the heliocentric distance to the Sgr dSph and its 
proper motion, and seven "noise" parameters. The initial 
angular position and heliocentric velocity of the Sgr dSph 
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are held fixed. Our priors for the disk scale length and scale 



height are based on Juric et al. (2008) while our priors for 



the Sun's peculiar velocity are based on values from |Binney | 



&; Merrifleld 



(1998). 



In Table 2 we present a summary of the statistics for 
the model parameters. In particular, we give the mean and 
variance of the marginalized PDF for each parameter. For 
convenience, we also list the prior used for each parameter. 



4.1 General properties of the Galaxy 

Table 3 presents a summary of the statistics for the ob- 
servables based on our MCMC analysis. For comparison, we 
include the predictions from the LJM model. We see that 
our analysis favors a somewhat lower value for the circu- 



lar speed than advocated by Bovy et al. ( 2009 ) though the 
discrepancy is only at the 1.5a- level. More importantly, our 
model is consistent with estimates for both the disk surface 
density and vertical force in the solar neighborhood in con- 
trast with the LJM model, which predicts values that are 
too high, implying that their disk mass is too large. 

In Figure [5] we show the LOSVD measurements from 
Tremaine et al.l (120021 and the inferred values from the 



model. Our model does a good job in matching the gen- 
eral value for the dispersion profile though the data show a 
clear peak at 200 pc and decline for larger radii whereas the 
model dispersion profile is relatively flat. The implication 
is that the simple Sersic-Prugniel & Simion model may be 
inadequate for the inner-most regions of the bulge. We see 
that the LJM model predicts LOSVD's that are 1.5-2 times 
higher than the measured velocities, which implies that their 
bulge mass is 2-4 times to high. 

The surface brightness toward the bulge is shown in Fig- 
ure [5] Both our model and the LJM model do an acceptable 
job in fitting the data. The inner rotation curve is shown 
in Figure [7] and again, an acceptable fit is obtained. Note 
that the LJM prediction rises more steeply at small \l\ as 
compared with our model and reaches a peak circular speed 
that is higher by a factor of 1.5. 

We show the outer rotation curve in Figure [8] In con- 
trast with the inner regions of the Galaxy, the LJM model 
and ours have a very similar shape. This point is further 
illustrated in Figure [9] where we show the full circular speed 
curve across a wide range in radii. For comparison, we also 
show the Xue et ah] ( |2008| rotation curve, though this data 
was not used in our analysis. The LJM model and ours 
yield remarkably similar circular speed curves even though 
the bulge-disk-halo decomposition of v c is very different. We 
note that the peak contribution to v c from the bulge is nearly 
twice that for our best-fit model, a result echoed in Figure 
[5] The peak contribution to v c from the disk is roughly the 
same in the LJM model as in our best-fit model, though in 
the LJM model, the peak occurs at R=10 kpc whereas in 
ours, it occurs at about 5-6 kpc. Recall that for an expo- 
nential disk, the peak in v c is at 2.2 disk scale lengths, so 
the implication is that the LJM disk has an effective scale 



length of 5 kpc, which is roughly a factor of two larger 
than the standard values for the Milky Way. 



4.2 The Sgr stream and halo triaxiality 

In Figure [lO] we show the angular position /3 and heliocen- 
tric velocity v r of stars that are presumed to be members 
of the Sgr stream. A few comments are in order. First, the 
M-giant data is considerably noisier than the SDSS fields 
though the latter only covers a portion of the leading arm of 
the Sgr stream. In Table 2, the stream-thickness (or noise) 
parameters are cmg = 6.9° and csdss = 3°. Thus, the anal- 
ysis accounts for the intrinsic scatter in the M-giant data by 
choosing a larger noise parameter. This choice has the effect 
of reducing the weight the individual M-giant data points as 
compared with the SDSS fields in the likelihood calculation. 

The PDFs for A p , B p , and are shown in Figure [TT] 
The inferred shape is still roughly oblate with the short axis 
nearly aligned with the Sun-GC line. In fact, the preferred 
model is more strongly aspherical; the best-fit values for 
A p and B p are even larger than were found in the analysis 
presented in Section 2. Note however that the uncertainties 
in these values is larger than before. 

Regardless of the model or the suite of constraints, the 
Sagittarius stream consistently points to a halo model that is 
oblate-triaxial with the short axis of the halo approximately 
coincident with the Sun-GC line. What properties of the Sgr 
stream drive the fit of the halo in this direction? In Figure 
|12| we show the isodensity contours for the total Milky Way 
model in the f3 — plane. Also shown are the principle 
axes of the halo projected on to this plane. Finally, we show 
positions of the M giant and SDSS fields along the stream. 
Roughly speaking, the stream (and presumably, the orbit 
of the Sgr dSph) are elongated parallel to the plane of the 
disk. In general, if isopotential surfaces in a particular plane 
are ellipses, then at least some of the particle orbits in that 
plane will be elongated perpendicular to the long axis of the 
potential. Thus, current observations of the Sgr stream may 
require a halo that is elongated perpendicular to the plane 
of the Galactic disk, as found in LMJ09, LM10 and in the 
present work. 



4.3 Mass distribution 

The Sgr stream inhabits the region of the Galaxy between 
20-60 kpc, which is roughly where the gravitational poten- 
tial transitions from disk-dominated to halo-dominated. The 
stream therefore has the potential to help break the disk- 
halo degeneracy and constrain the mass distribution beyond 
the disk. To explore this point, we re-analyze the model 
without the Sgr stream constraint. The result for the cu- 
mulative mass profile is shown in Figure [13] We see that 
uncertainties in the mass at r = 100 kpc are reduced from 
0.3 dex to 0.1 dex by including the Sgr stream,. 

As noted above, and illustrated in Figure [13] the mass 
distribution for our model and for the LJM model are very 
similar though the mass decomposition are very different. 
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Table 2. Statistical summary of results from the MCMC run described in Section 4. Column 1 
— parameter; column 2 — mean value; column 3 — variance; column 4 — type of prior used. 
Columns 5 and 6 give the lower and upper bounds for prior in the case of either a uniform or 
log prior and the width and mean for the case where the prior is a Gaussian. The units are 
km s -1 , kpc, Mq, Mq pc -3 , mas yr _1 for velocity, distance, mass, density, and proper motions 
respectively. 



Observation 


Observed 


LJM 


This Work 


v c (km s _1 ) 


244 ± 13 


220 


225 ±4 


Oort A ( km s _1 kpc -1 ) 


14 ±0.8 


13.0 


14.4 ±0.3 


Oort B ( km s" 1 kpc -1 ) 


-12.4 ±0.6 


-14.5 


-12.6 ±0.4 


E d ( M 0P c -2 ) 


49 ±9 


95 


49 ±4 


(27rG)- 1 |X,(l.lkpc)| (M 0P c -2 ) 


71 ±6 


100.1 


75 ±5 


M(r < 100 kpc) (10 11 M ) 


7 ±2.5 


7.96 


6.99 ±0.63 



Table 3. Mean values for observables from the LJM model and this work. 



In particular, the disk and bulge masses for our model are 
lower in our models by factor of 2-4. 



5 CONCLUSIONS 

We have performed a Bayesian analysis designed to model 
Milky Way that incorporates data for the Sgr stream to- 
gether with a broader suite of observational constraints. Our 



analysis allows us to not only infer the shape of the dark 
matter halo, but also its spherically- averaged mass distribu- 
tion. We find that our Bayesian approach, which uses single- 
particle orbits, yields results for the shape parameters of the 
halo that are consistent with the N-body-based approach of 
LM10. Moreover, the Sgr stream reduces the uncertainty in 
the Milky Way mass profile by a factor of ~ 2. 

Our analysis has not resolved the central question posed 
by LMJ09 and LM10: Why does the Sgr stream appear to 
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favour such a strongly triaxial halo? As discussed in the in- 
troduction, triaxial halo models favoured by LMJ09, LM10, 
and this work are outliers in the distribution of halo shapes 
predicted by the standard A-CDM cosmology. While the 
Galactic halo may indeed have this shape, we caution the 
reader that these analyses make use of a single stream. One 
should be suspicious that the halo model shares the same ap- 
proximate orientation as the orbital plane of the Sgr stream. 

In principle, the use of multiple stellar streams with dif- 
ferent orbital planes can improve the situation. Indeed, the 
GD-1, Orphan, and Monoceros streams have been already 
been extensively studied and, in some cases, used to con- 
strain Milky Way models (see, for example Koposov, Rix, 



fe Hogg (2010); Newberg et al. (2010); Penarrubia et al. 



(2005)). As discussed in Koposov, Rix, fe Hogg (2010), the 



GD-1 stream is a particular attractive feature since it is ex- 
tremely narrow and therefore has a simple internal structure 
as compared with the Sgr stream. However, the GD-1 and 
Orphan streams cover a significantly smaller angular extent 
than the Sgr stream and therefore may have limited value 
in constraining the potential. 

The tidal debris of satellite galaxies may, in fact, pro- 
vide other less direct clues about halo triaxiality. N-body 



simulations by Rojas-Nino et al. (2012) (see also |Penarrubia, 
Walker, & Gilmore (2009)) suggest that the morphology of 



tidal debris depends on the shape of the host galaxy's halo 
and in particular, the degree to which it departs from spher- 
ical symmetry. 

The results presented in this work and in LMJ09 and 
LM10 provide a tantilizing, if not unsettling picture of the 
Galactic halo. Clearly, further analysis and more data are 
required to establish what the shape of the dark halo is. 
Advances will likely come from a combination of improved 
modelling techniques, detailed N-body simulations and new 
data from the the next generation of surveys, such as Gaia 



Perryman| ( |2Q02| ) and LSST |Ivezic et al.|p008 ). 

We thank K. Johnston, R. Ibata, D. Hogg, and D. 
Foreman- Mackey for insightful comments and suggestions. 
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ences and Engineering Research Council of Canada through 
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Python programming language and the open-source mod- 
ules numpy and mat plot lib. 
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Figure 1. Marginalized PDFs for and 6 from the MCMC run described in Section 2.4 

(TP). Contours enclose 68%, 95%, and 99% of the probability in either the — or — 
planes (left and right panels, respectively). The black points indicate models from the MCMC run 
that lie outside the 99% contour. The stars show the best-fit parameters from LMJ09 (red) and 
LM10 (blue). The diagonal, horizontal, and vertical lines in the left-hand panel are for reference. 
Models along the diagonal line have = B§, and are axisymmetric with the symmetry axis 
coincident with the xt axis. Along the upper portion of this line the models are oblate. Models 
along the horizontal = 1 line are axisymmetric with the symmetry axis coincident with the 
z-axis (i.e., the symmetry axis of the disk.) 
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Figure 2. The PDFs for A p , B p and 6 for TD as described in Section 2.5 (upper panels) and 
TP (lower panels). For the upper panels, A<$> is transformed to A p using the expression A p = 
3 (A<£ — 1) + 1 and likewise for B. The contours, black points, and symbols are the same as in 
Figure 1. 
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Figure 3. The marginalized PDF in the R Q — Dq st from TD. The contours and black points are 
the same as Figure 1. The blue star shows the values assumed by LMJ09 and LM10. 
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Figure 4. The PDF in the fn — fii plane. The contours and black points are the same as Fig. 1. 
The black cross shows the results from the |Pryor et al.| (2007) analysis. The blue star is the value 
favoured by LM10. 
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Figure 5. Velocity dispersion profile as a function of projected radius. The thin blue line shows 
the velocity dispersion averaged over the chain. The shaded region gives the 68% credible region. 
The red dashed curve is the prediction for the LJM model. The black points are from |Tremaine| 
|et al.| ( |2002] ). 



© 0000 RAS, MNRAS 000, 000-000 



The Sagittarius Stream and Halo Triaxiality 15 




© 0000 RAS, MNRAS 000, 000-000 



16 N. Deg & L. Widrow 




Figure 7. Terminal velocity as a function of sinL The thin blue line and shaded region are the 
average and 68% credible region, respectively while the dashed red curve shows the prediction for 
the LMJ model. The black points are the data from Malhotra (1995). 
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Figure 8. Outer rotation curve as a function of projected radius R. The curves are as in Figure 5. 
The black points are from Brand & Blitz ( 1993 ) and the red points are from De mers &: Battinelli| 
( [2007] ). 
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Figure 9. Composite curcular speed curve as a function of radius. Lines and shaded regions 
show the average the 68% credible regions (upper panel) for the total circular speed curve (blue) 
and the contributions from the bulge (red), disk (green), and halo (magenta). The dashed curves 
(lower panel) are the prediction from the LJM model. Black points are from |Xue et al.| (2008). 
These points are not used in the fit but are shown as a consistency check. 
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Figure 10. Angular postion and line-of-sight velocity as a function of position along the stream. 
The upper panel shows the angular position perpendicular to the stream, /3, as a function of 
angular position along the stream, A. The lower panel shows the heliocentric radial velocity, v r , 
as a function of A. The thin blue line, shaded region, and dashed curve are as in Figure 5. The 
black points are the 2MASS M giants while the red stars are the |Belokurov et al.| (2006) SDSS 
fields. The red dashed line is the predictions for the best-fit model of LMJ09. While the LM10 
should be considered superior to the LMJ09 model, it produces the fit to the data using an N-body 
representation of the dwarf rather than the single particle approximation. Since we also use the 
single particle approximation, it is more appropriate to show that comparison. 
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Figure 12. Total Milky Way density and the Sgr stream in the /3 = plane. The green, blue, and 
red stars show the location of the Sun, GC, and Sgr dwarf respectively. The black points show the 
Sgr stream M Giants and the yellow points show the SDSS fields from Beloku rov et al.| (2006). 
The dashed curve shows the orbit of the dwarf in the Milky Way model with the parameters equal 
to the expectation values from Table 2. The solid blue, red, and black lines show the A, £?, and 
1 axes respectively. 
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Figure 13. Cumulative mass profile as a function of spherical radius r. Solid lines and shaded 
regions show the average values 68% credible region for the total (blue) and the separate con- 
tribution of the bulge (red), disk (green), and halo (magenta). The dashed red line shows the 
predictions of the LJM model. The upper panel shows M(r) with the full suite of constraints, 
include the data for the Sgr Stream. The lower panel shows M(r) when the stream data is not 
used. 
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